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We study the phase diagram and the commensurate-incommensurate transitions in a phase field model of 
a two-dimensional crystal lattice in the presence of an external pinning potential. The model allows for both 
elastic and plastic deformations and provides a continuum description of lattice systems, such as for adsorbed 
atomic layers or two-dimensional vortex lattices. Analytically, a mode expansion analysis is used to determine 
the ground states and the commensurate-incommensurate transitions in the model as a function of the strength 
of the pinning potential and the lattice mismatch parameter. Numerical minimization of the corresponding free 
energy shows good agreement with the analytical predictions and provides details on the topological defects 
in the transition region. We find that for small mismatch the transition is of first-order, and it remains so for 
the largest values of mismatch studied here. Our results are consistent with results of simulations for atomistic 
models of adsorbed overlay ers. 

PACS numbers: 64.60.Cn Order-disorder transformations; statistical mechanics of model systems, 64.70.Rh Commensurate- 
incommensurate transitions, 68.43.De Statistical mechanics of adsorbates, 05.40.-a Fluctuation phenomena, random pro- 
cesses, noise, and Brownian motion 



In nature there exist many modulated structures which posses two or more competing length scales. Such systems often 
exhibit commensurate-incommensurate (CI) phase transitions 01 Q], which are characterized by structural changes induced by 
competition between these different scales. Systems in this class that have received intense attention for many years include 
spin-density-waves ylHl in Cr and charge-density-wave systems [5] e.g. in "blue bronze". They are characterized by an order 
parameter (e.g., charge or spin density) that is modulated in space with a given wave vector q. In these systems the transition 
from a commensurate (C) state to an incommensurate (I) state is controlled by temperature and interaction with defects of the 
lattice 1 6]. Other systems of interest are materials which exhibit magnetic phase with spiral-like structures 0,0]. Finally, vortex 
lattices in superconducting films with pinning centers | 9] and weakly adsorbed monolayers 03 El on a substrate comprise a 2D 
realization of systems exhibiting CI transitions. In these systems the interparticle interactions are minimized by a configuration 
with a lattice constant a, while the substrate- adsorbed monolayer interaction is minimized by a configuration with lattice constant 
b usually incommensurate with a. 

The simplest theoretical model for a CI transition is the ID Frenkel-Kontorova (FK) model lfl2lfl3ll . It consists of a chain of 
particles interconnected by springs, representing an adsorbate layer, and placed in a periodic potential describing the effects of a 
substrate. The potential energy of a such system is given by 
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where a is the equilibrium lattice spacing of the chain in the absence of a potential and A is the stiffness. The potential function 
V (x n ) has periodicity b and can be approximated by a cosine function 

V(x n ) = V (l-cos(^-^y (2) 

When Vo is sufficiently small, the adsorbate lattice will be independent of the potential. This structure is called a "floating" 
phase and the lattice spacing a = lim^oo (x n — xo)/n of the adsorbate lattice can be an arbitrary multiple of the substrate 
periodicity. In general, the floating phase is incommensurate for almost all values of the ratio a/b. 

In the opposite limit, when the potential is very strong, one can expect the lattice to relax into a commensurate structure where 
the average lattice spacing of the adsorbed atoms is a simple rational fraction of the period b. In the I phase, close to the CI 
transition, it is energetically more favorable for the system to form C regions separated by domain walls in which the springs 
are stretched or compressed and the commensurate registry with the potential is lost. These domain walls are usually called 
discommensurations and the corresponding region in the phase diagram within the I phase will be referred to as a modulated 
(M) phase. A positive (negative) discommensuration leads to a reduction (increase) in the density of adsorbate atoms and these 
regions are referred to as light (heavy) walls. 

The CI transition can be determined by examining the behavior of the winding number Q = a/b as a function of Q = a/b 
for fixed Vo. In the continuum version of the FK model the CI transition is second order with a correlation length, identified as 
the domain wall separation, which diverges logarithmical near the transition. If Vo is larger then a certain critical value, V c , the 
system will be commensurate for all the values of O. The function Q(Q) has a staircase type of appearance and thus for Vo < V c 
it is called an incomplete devil's staircase and for Vo > V c a complete devil's staircase lEll4ll. If there are discontinuities or first 
order jumps between commensurate states it is called a harmless staircase lfl5ll . 

The FK model can be extended to two dimensions to describe, for example, adsorbed layers on crystal surfaces. In its simplest 
version the 2D FK model (3l describes the adsorbate interactions by a pure harmonic potential in a periodic pinning potential. 
Although the FK model takes into account topological defects in the form of domain walls it leaves out plastic deformations of 
the layer due to topological defects such as dislocations. These defects are particularly important when the CI transition occurs 
between two different crystal structures or in presence of temperature fluctuations or quenched disorder. 

They are automatically included in a full microscopic model involving interacting atoms in the presence of a substrate po- 
tential. However, the full complexities of the microscopic model limits the actual numerical computation to systems of relative 
small sizes. The size effect could be very strong for CI transition involving extended domains or topological defects. Recently a 
phase field crystal model was introduced lialljl that allows for both elastic and plastic deformations in the solid phase. In this 
formulation a free energy functional is introduced which depends on a density averaged over microscopic times scales, <3>(r, t). 
The free energy is minimized when is spatially periodic (i.e., crystalline) in the solid phase and constant in the liquid phase. 
By incorporating phenomena on atomic length scales the approach naturally includes elastic and plastic deformations, multiple 
crystal orientations and anisotropic structures in a manner similar to other microscopic approaches such as molecular dynamics. 
However, the PFC method describes the density on a diffusive and not the real microscopic times scales. It is therefore com- 
putationally much more efficient. Thus this model should provide a suitable description of the CI transition when topological 
defects, domain walls and dislocations are present. 

In this work we extend the 2D phase field crystal (PFC) model 0[n|] to include the influence of an external pinning potential. 
Such a model should provide a suitable continuum description of lattice systems such as weakly adsorbed atomic overlayers or 
2D vortex lattices with pinning. The pinning potential is chosen such that it induces CI transitions between ground states of 
different symmetries in the model. The outline of the paper is as follows. We first define the model and carry out an analytic 
mode expansion to determine the crystal structure and the location of the CI transition lines. We analyze these transitions as 
a function of the lattice mismatch and strength of the potential. Following this we carry out a numerical minimization of the 
full free energy functional which gives good agreement with the theoretical predictions and provides details on the nature of 
the topological defects near the transition. Finally we discuss the nature of the CI transitions and the relation of our results to 
atomistic simulations for adsorbed overlayers. lfl8ll. 

II. THE PHASE FIELD CRYSTAL MODEL 

In the phase field crystal model 00, the free energy functional is written as 

r = f*(^L* + u * + *G&)*), (3) 

where G(V 2 ) = A {q$ + V 2 ) , and its eigenvalues can be related to the experimental structure factor of e.g. Ar E3. Eq. © 
describes a crystal with a lattice constant of 2n/q , while the elastic properties can be adjusted by A, u and qo. For numerical 
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calculation it is convenient to rewrite the free energy in dimensionless units as 
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In these units the free energy becomes 
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(l + V 2 ) 2 . Since ip is a conserved field it satisfies the following equation of motion, 
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where ( has zero mean, (( (xi, ti) C (^2? ^2)) = VV 2 5 (x\ — X2) S (ti — and P is a constant. Equations Q and (|6|) 
have been used to study a variety of phenomena involving elastic and plastic deformation including grain boundary energies 
between misaligned grains, buckling and dislocation nucleation in liquid phase epitaxial growth, the reverse Hall-Petch effect in 
nano-crystalline materials, grain growth and ductile fracture 11711. 

In 2D the free energy in Eq. in the absence of external potential is minimized by three distinct solutions for the dimension- 
less field ip\ constant, stripes and dots. For the purposes of this work only the latter solution is relevant. This solution consists of 
a triangular distribution of density maxima corresponding to a crystalline phase. In general this solution can be written as 



^0) = (ir] 
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where G nm = nb\ + m&2, &i and 62 are the reciprocal lattice vectors and ip is the average value of ip. For a triangular lattice b\ 
and 62 can be written as 
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where a t is the distance between nearest-neighbor (local) maxima of ip (corresponding to the "atomic" positions). The ampli- 
tudes a m?n and lattice spacing a t are determined by minimizing the free energy functional. For simplicity it is useful to first 
consider a simple one-mode approximation (OMA) in which only pairs (n, m) that correspond to |G n?m | = 2n/ (a t V3/2) are 
retained. In this limit ip can be written, 

ipt = A t [cos(q t x) cos(q t y/Vs) - cos(2q t y / Vs) / 2] + ip, (9) 
where A t is an unknown constant and q t = 2k /at. Substituting Eq. (|9|) in Eq. © and minimizing with respect to A t and q t 

gives, A t = 4(^ + (-15r - 36^ 2 ) 1/2 / 3 )/ 5 > g t = >/3/2 and 
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As shown in Ref.[n|this approximation is valid in the limit of small r. In the next section an external pinning potential with 
square symmetry and incommensurate with will be introduced in the PFC model. As will be shown , depending on the strength 
of the external potential and the lattice mismatch, a commensurate incommensurate transition would occur. 

III. PHASE DIAGRAM WITH EXTERNAL POTENTIAL 

In the PFC model , an external potential V is introduced by adding a term coupling V linear to ip in the free energy functional 
given in Eq. i.e., 
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In this study, we consider an external potential of square symmetry which is distinct from the symmetry of the triangular 
lattice in the absence of the external potential (e.g., Eq. |9j). 

V = Vo cos(q s x) cos(q s y). (12) 
where q s = 2ir/(a s a/2). We define the relative mismatch S m between the external potential and adsorbed monolayer as 

5 m = l- 2ir/a s (13) 
With the above external potential, the PFC model could describe for example an adsorbed layer on the (100) face of an fee crystal 

To understand the influence of V on the minimum energy solution, we can again Fourier analyze the equilibrium density. 
Taking into account both the intrinsic triangular symmetry and the external potential of square symmetry, the system is best de- 
scribed by a combination of hexagonal and square modes. In this case the corresponding hexagonal- square mode approximation 
(HSMA) can be written as 



ip hs = Ai cos(q s x) cos(q s y) + A 2 (cos(2q s x) + cos(2q s y)) 

+ A t [cos(q t x) cos(q t y/V3) - cos(2q t y/V3)/2]+^. (14) 

This Fourier expansion includes the basic mode for the triangular lattice as done in Eq. |9|and the first two harmonics for the 
commensurate modes with square symmetry. In general, they describe an incommensurate phase corresponding to an triangular 
phase distorted by the square symmetry external potential. However, when the mismatch is sufficiently large and/or the strength 
of the external potential is sufficiently strong, the solution that minimize the free energy will correspond to a vanishing value 
for A t and a commensurate phase. 

The free-energy density in the HSMA approximation is given by 

F '"' s - h A * + h A > + * (I + 's + f - b- + b-) ^ - ^ 

+ |-4? (I A? + 4) + A\ {^Al + \ A2 fj + Ml + F'/S. (15) 
The coefficients A t , A\ and A 2 are unknown and must be chosen to minimize the free-energy density. 

The analytic expressions for the free energies can now be used to obtain the phase diagram of the pinned PFC model as a 
function of the pinning strength Vq and the mismatch S m . To this end Eq. JT5l can be used to determine the critical value 
of Vo,V c , at which the CI transition occurs, identified by the point where the amplitude of the triangular phase vanishes, i.e., 
when At(Vo,5 m ) = 0. The results are shown in Fig. da). Within the HSMA the CI transitions is discontinuous and becomes 
continuous only for infinite S m . In the next section we will compare the analytically obtained phase diagram in the HSMA 
approximation with a full numerical minimization of the total free energy. 



IV. NUMERICAL RESULTS 
A. Energy Minimization 

While the HSMA yields a good qualitative understanding for the phase diagram describing the CI transition, it is not quan- 
titatively accrate, particularly near the transition, since the HSMA cannot describe the structure of the domain walls and other 
possible topological defects. In this section, we desribe the full numerical investigation of the CI transition in the PFC model. 
This is obtained by direct minimization of the free energy functional without any assumed form for the density field. 

In the presence of the external potential, the equation of motion for ip becomes, 



^=V 2 (co (V 2 ) ^ + V 3 + V) + (16) 

The minimal energy numerical solutions for the equilibrium states of ip were obtained by solving SF/Sip = using a simple 
relaxational method similar to the usual molecular dynamics annealing scheme. The noise term ( in Eq. fi"6l is only used as an 
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FIG. 1 : (a) The phase diagram calculated analytically using the HSMA approximation for the free energy. The region denoted by C is the 
commensurate phase and I denotes the incommensurate phase, (b) The corresponding phase diagram determined numerically as explained in 
Sec. |IV| The dashed line corresponds to a crossover regime between the fully incommensurate I and the modulated M phases. See text for 
details. 



annealing tool to escape from any metastable states. It is applied for a limited period of time only and then set back to zero. For 
the remainder of the time, Eq. (fT6fc was solved using a simple Euler algorithm, i.e., 

^n+i,i,j =V^n,i,i + ArV 2 ( (r+ (l + V 2 ) 2 ) tl>n,ij + 4>n,i,j + V n ,i,j ) , (17) 

where the Laplacian operator V 2 represents the lattice second order derivatives. A so-called spherical Laplacian approximation 
fl7ll was used for V 2 . For a thorough discussion on solving solving differential equations numerically using stencils, see 
Ref. |l9]. 

Eq. fTTt was solved on a 512 x 512 grid with the spatial discretization dx = 1 and time step At = 0.02 using periodic 
boundary conditions. The parameters r and %[) were chosen to correspond to a crystalline region of phase space, i.e., r = —1/4 
and -0 = —1/4. However, a hexagonal lattice cannot be fitted in a square geometry and in order to satisfy the periodic boundary 
conditions, the lattice will be distorted and it may exhibit domains of different orientations separated by walls, thereby giving 
a free-energy density higher than that of the ground state. The value of dx also influences the free-energy density when lattice 
constant is small. We have checked the finite-size effects for our lattice in the absence of pinning potential, the relative difference 
between free-energy density F n / S and the OMA is only about 0.56%. 



B. Phase Diagram and Ground State Configurations 

The free energy density and the structure factor *S'( | fc| ) = S(k) = |^(|/c|)| 2 were calculated numerically for several different 
values of 2ir/a s as allowed by the periodic boundary conditions. In Fig. |2|a set of configurations for the model with increasing 
amplitude of the pinning potential Vq are shown. As expected for V p = the ground state has perfect triangular symmetry. With 
increasing amplitude of the square pinning potential the configurations become spatially distorted and eventually the system 
undergoes a transition to a square lattice. 

We define the position of the CI transition by analyzing the structure factor S(k). The transition occurs when the peaks of the 
structure factor that correspond to the square lattice begin to split as shown in Figs[3]and|4l 

Next the approximate analytic solutions for the free energy are compared to the numerical ones as a function of pinning 
potential for fixed mismatch. As seen in Fig.|5]the free-energy density decreases as a function of the pinning strength. For small 
values of S m the transition to square symmetry occurs at relatively weak pinning strengths and is first-order. For larger values 
of 5 m the free energy decreases slowly and the CI transition appears continuous (see Fig. 0a)). However, we have checked 
the finite-size dependence up to systems of linear size 1024 and find that for up to at least 5 m = 0.52 the transition remains 
first-order. We note that a first-order CI transition is in fact consistent with theoretical predictions based on models of interacting 
domain walls on a pinning potential with square symmetry [20]. In Fig. |5]we can also see that the analytical HSMA gives 
reasonable agreement with the numerical results, especially for large misfit where domain structures are not as predominant. 
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(a) V = 0.00 (b) Vb = 0.02 (c) V = 0.04 (d) V = 0.05 



FIG. 2: Snapshots of ground states and the corresponding structure factors (lower figures) showing the transition of the system from a triangular 
lattice for zero pinning to a square lattice for S m = 0.14 (only a small part of the lattice is shown here). 




(a) Vb = 0.04440 (b) V = 0.04455 



FIG. 3: Splitting of the peaks of the structure factor with decreasing pinning potential in the vicinity of the CI transition for 5 m = 0.14. The 
splitting is discontinuous, suggesting a first-order like transition. 

Finally, we determined the onset of the modulated (M) phase within the I phase by the value of Vq where the peaks in the 
structure factor corresponding to the triangular symmetry split in multiple peaks. In the present model the M phase corresponds 
to regions commensurate with the pinning potential, but separated by heavy domains walls of excess of local density (see also 
Fig. The dashed lines indicates the threshold value where domain walls appear in the system, shown in Figs. [l]and|5ja)-(c) 
(in the insets) by dashed lines. 
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(a) Vo = 0.4275 (b) V = 0.4455 



FIG. 4: Splitting of the peaks of the structure factor with decreasing pinning potential in the vicinity of the CI transition for 5 m = 0.52. The 
splitting is much smaller than in Fig. 3 but remains discontinuous. See text for details. 




V o 

(c) S m = -0.01 



FIG. 5: Comparison of the free-energy densities between analytical (HSMA) and numerical results for different values of <5 m . The inset 
represents intensity of the peaks of the structure factor corresponding to the reciprocal vectors of magnitude k = 2it /a s . The continuous 
vertical lines represent the position of the CI transition, while the dashed lines marks the cross-over regime in the I phase. 



C. Voronoi Analysis of Domain Walls 

It is of further interest to analyze the nature of the spatial configurations and domain walls in the I phase and in the transition 
region. A convenient way to quantify the defects is to use the Voronoi analysis 12 ll 12211 . By definition, the Voronoi cell of any 
particle contains all the points that are closer to it than any other particle, i.e. it defines the Wigner-Seitz cell for each particle. 
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The resulting cells are polygons with N sides that represents the number of nearest neighbors (NNs). In the present case this 
is particularly useful since very close to the transition region there is some variation in the positions and sizes of the density 
maxima, identified as effective "particles" in our phase field crystal model. In Fig. [6] we show the results of Voronoi analysis of 
our data as a function the pinning strength for a fixed mismatch S m = 0.14. For weak pinning, Voronoi cells with N=6 dominate 
as expected. Closer to the CI transition line (shown with a vertical line) contributions from N=5 becomes important and in the 
immediate vicinity of the CI transition virtually all Voronoi cells have N=4. The Voronoi cell analysis can be used to identify 
the location and nature of the domain walls. Close to the CI transition, the system is composed of regions commensurate with 
the pinning potential but separated by heavy walls. These configurations appear spontaneously and have been observed also 
in Monte Carlo simulations of overlay ers of atoms interacting via the Lennard- Jones potential adsorbed on the (100) face of 
an fee crystal 111 811 . In Fig. [7] we show results for the present model in the corresponding region. The density of maxima is 
smoothed with a gaussian function with width of approximately 3a s . As can be seen in Fig. |7(a)| the system indeed comprises 
commensurate regions separated by walls. The regions of high density indicate "heavy" domain walls with excess of "particles" 
compared to the commensurate regions . In these regions the "particles" are plotted with black circles while in the commensurate 
regions with white circles. In Fig. |7(b)| we show the contribution of different Voronoi polygons in this state. The fraction of 
polygons with N=5 (ps) is 0.53, while the fraction with N=4 (P4) and N=6 ($>$) are 0.16 and 0.29, respectively, which agrees 
well with results in Ref.fl8l 
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(a) 



(b) 



FIG. 7: (a) Example of modulated phase obtained from system characterized by 5 m = 0.14, Vb = 0.043. The dark areas represent the I 
regions, i.e. heavy walls, while the light areas are C regions. The positions of the 'particles' are superimposed on the density maxima, (b) The 
contribution of Voronoi polygons in this state. See text for details. 
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V. DISCUSSION AND CONCLUSIONS 

In this work we have considered the recently developed crystal phase field model fnll in the presence of an external pinning 
potential. As the model naturally incorporates both elastic and plastic deformations, it provides a continuum description of lattice 
systems such as adsorbed atomic layers or 2D vortex lattices, while still retaining the discrete lattice symmetry of the solid phase. 
The main advantage of the model as compared to traditional approaches is that despite retaining spatial resolution on an atomic 
scale its temporal evolution naturally follows diffusive time scales. Thus the numerical simulations studies of the dynamics 
of the systems such as approach to equilibrium can be achieved over realistic time scales many order of magnitudes over the 
microscopic atomic models. We have taken advantage of this and considered the full phase diagram of the model as a function 
of the lattice mismatch and pinning strength, both analytically and numerically. A systematic mode expansion analysis has been 
used to determine the ground states and the commensurate-incommensurate transitions in the model. Numerical minimization 
of the corresponding free energy shows good agreement with the analytical predictions and provides details on the topological 
defects in the transition region. In particular, we find that the transition remains discontinuous for all values of the mismatch 
studied here. We have also performed a detailed Voronoi analysis of the domain walls throughout the transition region. Our 
results are consistent with simulations for atomistic models of pinned overlay ers on surfaces. 

A particularly interesting application of the present model is to pinned systems which are driven by external force. Examples 
of such systems include driven adsorbed monolayers ll23ll . driven charge density waves ll24ll and driven flux lattices ll25ll . Work 
in these problems is already in progress. 
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